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ABSTRACT 

The quality and ef f ectiveneea of decision inakina relates 
directly to the quality of data , analytical tools available and 
sophistication of mathematical model adopted by decision-makers. 
A close descT-lption of Earth’s relief is needed for several 
purposes by administrators, planners and Resource scientists. 
Such a requirement cannot be met in normal working, except by 
way of supplementation through various means. Digital Terrain 
Model (DTM) defined as registration of array of coordinates Ja 
one such means. The accuracy of DTM depends on density versus 
sophistical ioti of mathematical model Z » F(x,y) approximating the 
surface . In this investigation, an attempt has been made to 
evolve a unique system of DTM for volumetric computations. 


In Civil Engineering Projects Volumetric computations play a 
unique role. Particularly in highway design phase, such 
determinations are essential for both alignment and cost 
estimation. Volumetric determination through DTM offers an 
advantage to designer to compare several alignments without 
recourse to further measurements. Apart from the academic value, 
this work will have great use in practical methods for 
planning and execution of Civil Engineering projects. 
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CHAPTER I 


INTRODUCTION 


1.1 Introduction 


The word Photofirannetry is derived from Greek words 
'Photos', 'erasuna' and 'metron' , meaning 'light', 'that which is 
drawn or written’ and 'to measure* respectively. The word 
'photogrammetry ’ literally means measuring from photographs. 

The last few years have seen a change in photogrammetry from 
graphic products to digital, brought out by ubiquitous computers. 
Another change is on the horizon now — Digital Terrain Nodelling 
(DTM), which need structured data. "Digital Terrain Modelling is 
self explanatory -- a DTK is the digital representation of the 
terrain i.e., a representation of ground surface in terms of 
X,Y,Z coordinates. X and Y define planimetric position of a 
terrain point, Z its height. The various uses of DTK include 
(1) Production and maintenance of cartographic data base (2) 
Orthophoto production (Orthophoto is photograph showing images of 
objects in their true or orthographic positions) (3) Engineering 
planning such as highways, railroads and power lines etc., (4) 
Simulation (5) Keteorology and Navigation. 

The data can be organized as a matrix array of coordinate 
triplets or as equations of surface defined by polynomials or 
Fourier series. Characteristics other than elevation such as 
land value, soil type, depth to bedrock and land use etc., may 
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also be included in DTM. Uhen only elevation is considered it is 
normally known as Digital Elevation Hodels. 

1.2 Methods of Surface Representation 

The data for digital terrain model can be stored as an array 
shown in Fig. 1.1. The grid spacing in X and Y direction remains 
constant and only elevation (Z-coordinate) of the data points are 
recorded. Knowing the grid spacing and the number of rows and 
columns in the grid exact X and Y coordinates of each elevation 
point can be easily derived. However, since the placement of 
data points is completely independent of surface features, this 
method of data distribution may result in the omission of 

critical surface features. The accuracy of surface 

representation will depend largely on grid size and surface 
roughness. 

The second method of representation of terrain is in the 
form of a polynomial. Earth’s surface may be represented by 3, 
4,8,12 or 16 parameter bicubic polynomial in X,Y coordinates. A 
16-parameter polynomial could be 

Z = A + BY + CY^ + DY^ + EX + FXY + GXY^ + HXY^ + IX^ + JX^Y + 

KX^Y^ + LX^Y^ + MX^ + NX^Y + PX^Y^ + QX^Y^ 

where A,B,...,P,Q are 16 parameters. The accuracy of 
surface depends on density of data points collected, 
sophistication of the mathematical model Z = ffX.Y) and type of 
terrain. This is shown fig 1.2. 

Truly speaking for complex surfaces the continuous function 




piCsI-i SURFACE REPRESENTATION BY POLYNOMIAL 
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Fi'gl3 SURFACE REPRESENTATION BY FOURIER SERIES. 
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governing the aurface will be a polynomial of the type 

Z » Jaj^xV. 

Surface rapresentat ion using Fourier series is a very 
advanced method of approximating the surface. The terminology 
associated with Fourier series is derived largely from Electrical 
Engineering and time series analysis. 

In the simplest example, all harmonics in one direction have 
zero amplitude, and only one amplitude in other direction has an 
amplitude greater than zero. The resulting surface resembles of 
a sheet of corrugated iron roofing or a sheet of parallel 
ripples, 'or a still body of water as shown in Fig. 1.3. 

The mathematical background of Fourier aeries is excellently 
given by John C. Davis [6]. 

1.3 Need and Objective of Uork 

Earth’s relief is needed for diverse purposes as in 
contouring on maps, determination of intervisibility of points, 
terrain simulation, slope analysis. Earth work calculations and 
various Engineering applications. 

Computation of volumes for highways from conventional 
methods is tedious and expensive. Alternative alignment 
selection tlirough ground survey measurements is expensive and 
time consuming process. A rectangular DTH generated by any 
technique is useful for the above objectives. 

The objective of the thesis is to describe 

(1) various methods of data collection 

(2) to explain the method to generate DTM for flat terrain 
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and undulated terrain . 

(3) to diacuss the method of progreaaive aamplinfi to generate 
cartographic data base 

(4) to develop a computer program for highway volume calculation, 

given a DIM and design parameters of highway. 

1.4 Organization of the Thesis 

Previous works on DTM are discussed in Chapter II. Chapter 
III describes the various sources of data collection for DTM. 
Chapter IV is devoted to mathematical formulation and Chapter V 
describes volume calculation methods. Results are summarized in 
Chapter VI. Chapter VII sums up the conclusions drawn from the 
present study and offers certain suggestions for further study. 



CHAPTER II 


LITERATURE REVIEU 


This Chapter is devoted to a survey of existing literature 
on DTM . A review of the existing DTIfs outside India and the work 
done in India in this realm are presented here. The literature 
available on the application of DTM to volumetric computations is 
very limited. But, in this thesis the literature available on 
datji acquisition, interpolation methods and applications are 
interrelated to evolve a unique system of DTM for volumetric 
computations. 

2.1 Introduction 

DTM was the brain child of Charles L. Miller at Massachusetts 
Institute of Technology ca 1955. The objective was to expedite 
highway design by digital computation based upon photogra- 
mmetrically acquired terrain data. Instead of measuring levels 
along cross sections of a given alignment, X,Y,Z coordinates of 
points in the band of interest are measured ("band of interest is 
a strip of ground in which it is known that the final horizontal 
alignment of road will pass"). The design of highway can take 
place with the aid of computer. The parameters for any alignment 
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may be fed in to the computer, volumes, mass movements etc., can 
be computed by each version which enables the desifiner to compare 
several alifinments without recour-iie to further measurements. 

2.2 Brief Description of some existing DTMs 


2.2.1 French DTM — "semis des points” 


Sampling pattern 


: morphologic points and lines 


Interpolation methods : moving average - to calculate the 

height Z of a point with known 
X,Y, a square is defined around 
this point. A polynomial of 
second degree is fitted using 
least squares. The function 
value of the computed parboiled 
at (X, Y) is taken as the height 
of this point. 

Application : highway design 

2.2.2. Finnish DTM 


Sampling pattern 
Interpolation method 


regular triangular grid, 
linear interpolation. 


Application 


highway design 
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2.2.3. Hertfordshire Counl.y Council 0TM 
Sampling pattern : Contours 

Interpolation : linear interpolation 

Application : highway design 


2.2.4. Czechoslovakian DIM 

Sampling pattern ; depends on the type of terrain 

Interpolation : Moving averages 

Appl ical.ion : highway design 

The work of Gogia, Survey of India, Chandigarh [5] reveals 
that for flat terrain least squares technique can be used to 
generate DTM. He suggested that bilinear interpolation can ' be 
adopted for undulated/hilly terrains. 

NaitanJ, in his work, "Photogrammetr ic methods for Rail 
route surveys and Profiling” [10] used 1:10000 maps and second 
order photogrammetric equipments, viz., wild B-8 for mapping. He 
suggested for corridor area study DTM is much useful. But this 
aspect was not discussed in his paper. 

Makarvic, in his paper, "Progressive sampling for digital 
terrain models" [8 ] stated that the method of progressive 
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sampling of discrete points arranged in a regular grid of varying 
intensity provides a powerful tool for acquisition of digital 
terrain model data. 

Ghazal 1 1 t 4 ] has done extensive investigation of the 
performance of progressive aampllng. He suggested numerical 
values of threshold based on nature of the terrain. 

Volumetric computation may be carried out using one of the 
two methods. 

(1) Crosssection method 

(2) The unit area or borrow pit 


The Crosssection method 

When Crosssection are used ground profiles are taken 10-50 
meters normal to present center line. The formula that is 
commonly used is end area formula . 

In this end area formula, volume is assumed to be equal to 
average area of the two ends multiplied by the distance between 
them. 

V = (Ai + A2) L/2 

where Aj^ and A 2 sre areas of two end bases and 1. is the 
distance between them. Because of its simplicity this formula is 
almost universally used in computing earthwork. 

This is illustrated in fig 2.1. 
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Fig 2.1 AVERAGE END AREA METHOD 
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Borrow PIT method 


Uhen 

the 

work 

covers a broad area 

rather 

than a narrow 

strip, as 

is 

the 

case in grading for 

a 

large 

building 

or 

industrial 

yard , or 

when a borrow pit is 

not adjacent to 

the 

route survey, 

it is 

frequently desirable 

to 

lay 

out 

a system 

of 

rectangles 

(usually 

squares) say 50 ft 

or 

100 

ft 

on a side 

depending 

on 

the roughness and topography 

of 

the 

ground 

and 


precision with which quantities are required. 

The volume of a square or a rectanfiular truncated prism can 
be calculated by 

V = A C hj^ + h2 + hg + h4)/4 

A is the area of the right section 

hi _4 are cuts at the four corners of the prism. 

Similarly, volume of the truncated triangular prism is 
calculated using the formula 

V = A(hi + h2 + h3)/3 

This method is illustrated in Fig. 2. 3. 

2.3. Proposed Method 

From the available literature, to achieve the objectives of 
the thesis the following method is proposed. 

(1) Data is collected from maps 

(2) Using weighted average method, the data is airanged into 

grids 



16 


(3) For further densif ication, least square interpolation 
is used for flat terrain and bilinear interpolation is 
used for undulated terrain . 

(4) Progressive sampling is done to locate new points 

(5) For volume calculations cross section method is used. 
This is schematically shown in fig2.3. 

A serion of computer programs coded in Fortran IV and 
implemented on DECS7STEII 1090, based on the above approach, makes 
it possible to achieve the objectives of the thesis. 
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CHAPTER III 

DATA ACQUISITION 

The various methods of data acquisition for digital terrain 
models are described in this Chapter. 

3.1 Sources of Data Collection 

There are three possibilities to obtain X,Y,Z coordinates of 
terrain points. 

(a) Ground Survey 

Cb) Photogrammetric Nodels 

(c) Topographic Maps 

3.1.1 Ground Survey 

In ground survey tacheometry or leveling is adopted to 
collect the data. Its a tedious and time consuming task. 

3 . 1 .2 Photogrammetric Models 

These are commonly used to collect data. The data can be 
collected using stereo plotters. Stereoscopic plotting 
instruments (commonly called stereo plotters or simply plotters) 
are instruments designed to provide rigorously accurate analog 
solutions for object positions from their corresponding image 
positions on overlapping pairs of photos. In aerial photography, 
the photographs are taken in such a way that there exists 60% 
overlap in longitudinal direction and 20% in transverse direction 
as shown in Fig. 3.1. To form a three-dimensional model of the 
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LATERAL 
OVERLAP 20V. 


Fig3.1 OVERLAPPING PHOTOGRAPHS- 



Fig 3.2 INTERIOR ORIENTATION. 
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Fig 3.3 STEREOMOOEL FORMATION . 

Page 




?2 

terrain an overlapping pair of diapoaitivea or tranaparen-ci ea 

are placed into two projectora of photogrammetr ic nachinea and 

1 

are reprojected. Thia procedure ia called interior orientation 
and ia ahovn in Fig. 3. 2. 

In the proceaa of reprojection the two adjoining projectora 
ahould be provided with auch movementa aa to correapond with the 
poaaible camera poaitiona at the time of photography ao that by 
their mutual adjuatmenta, the aame relative poaition of adjoining 
photographa aa during photography can be obtained. The proceaa 
ia called relative orientation and createa, in miniature a true 
three-dimenaional atereo model of the overlap area. 

After relative orientation, abaolute orientation is per- 
formed. In this proceos the model ia brought to desired scale 
and leveled with respect to reference datum. Fig. 3. 3 illustrates 
a stereo model created by interior, relative and abaolute 
orientation of a atereo pair in the projectors or a Balpex 
plotter. For measurements, the poaition of any point is 
determined by bringing a reference mark into contact with model 
point. Planimetric positions (X,Y) of points are plotted by 
means of a pencil located vertically beneath the reference mark, 
and elevation CZ) are read directly from dial and which records 
up and down motion of platen, the dial having been indexed to 
ground control during abaolute orientation. Fig. 3. 4 describes the 
principal components of a typical direct optical projection’ 
stereo plotter. 



8 



Fig 3.4 PRINCIPAL COMPONENTS OF STEREOPLOTTER 
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PRINCIPAL COMPONENTS OF STEREOPLOTTER 


(1) Main frane 

(2) Reference table 

(3) Traclnfi table 

(4) Platen 

(5) Guide rode 

(6) Projectors 

(7) Illuninatlon lamps 

(8) Diaposltlves 

(9) Leveling screws 

(10) Projector bar 
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a .l.at'a) Data ColXeetioft fcom Stereo Plotters 

Uith the analytical plotter two photographs are 
at ercoscopically viewed with a special optical system (Fig. 3. 5). 
In the centre of the stereo photo there is a measuring mark. The 
coor-dinates of terrain point, where the measuring mark can be 
seen are stored in the computer as model coordinates. If the 
operator changes the coordinates, normally by turning hand 
wheels, the computer moves the photos in front of the optical 
system so that are can see the measuring mark at the point with 
the new coordinates. The task of the computer is, besides moving 
photos, putting in the paramelors , doing the interior, relative 
and absolute orientation and error correction. The task of the 
operator is to produce input data for the different processing 
models. According to the mode, the the operator positions the 
measuring mark on single point or follows lines. The 
registration of the coordinates will be done by the operator or 
automatically in fixed time or space intervals. The control of 
the servos and thus the motion of the photos is automatically 
carried out by computer. 

Uithin each stereo model, an X-Y horizontal coordinate 
system is established with X-axis in the general direction of 
proposed highway alignment. Profiles in Y-direction are taken at 
regular intervals, and the Z-elevation along these profiles is 
recorded either at regular intervals of Y or at significant 
breaks in terrain slope, or both. 




Fig 3.5 PRINCIPLE OF ANALYTICAL SYSTEMS 


0 M PU T E R 
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3.1.3 Topofiraphic Hapa 

A nap is a tuo-dinensional representation of an area. Haps 
are inportant to photo^rannetrist for they are compact and 
efficient means of expressinfi relationships and details. The two 
ways of data collection from a map are by diftitizins contours 
(this will be discussed in sampling patterns) and from spot 
heights available from the map. The pattern of points on a map 
may be classified into (a) points regularly spaced on a grid or 
network (b) points scattered at random (c) points grouped in 
clusters. This is illustrated in Fig. 3. 6. 


For 

the present study, the spot 

heights 

and contours 

from 

Map (restricted) are 

collected . 

The point distribution 

is 

scattered 

at random. 

The scale of 

the map 

is .1:10000. 

Spot 

he! ghts 

are collected 

taking left 

bottom 

corner as origin. 

Heights 

are collected 

in meters. The 

matrix 

of data is used 

for 


the analysis of undulated terrain. For flat terrain analysis, 
the data is collected from reference [9 ], as the real world data 
was not readily accessible. 

In general one or combinations of these sources is used. 
The problem of which points on the ground surface should be 
selected to represent the terrain can be split up into the choice 
of a pattern at point registration (sampling pattern) and an 
adequate sampling density. 
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Fig 3,6 PATTERN OF POINTS ON MAPS. 








3.2 Pattarna 


There are numeroue possible sampling patterns, the most 
frequently used ones are discussed briefly here. 

a) Homogeneous Grids 

A rectangular or square grid is used to define the position 
of every point where a height has to be measured. The obtained Z 
coordinates are stored in matrix form. 

Advantages : The photogrammetric model can effectively be 
used for the sampling process. The sampling does not require a 
great deal of expertise, and can be automated. A simple 
interpolation and retrieval system is possible. A DTM system 
based on a regular grid can be implemented on a relatively small 
computer . 

Disadvantages : The cost-accuracy relationship is not likely 
to be good. (In uniformly shaped areas too may points will be 
measured and in strongly varying terrain too few, in general). 

b) Contour Lines 

Points are sampled along lines with a constant height, 
either in a photogrammetric model or on a topographic map. The 
obtained X and Y coordinates of the points are stored sequin- 
tali for each given height value. 

Advantages : This pattern of registration is much more 
related to the terrain features than a homogeneous grid. There- 
fore they have the potential to achieve a good accuracy economy 
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relationship. Digitizing contours or plotting contours in a 
photogranunetric instrument does not require a high degree of 
expertise. 

Disadvantages : Continuous sampling is less accurate than 
point by point sampling. The sampling process can hardly be 
atomized. Data compression procedures will be necessary, in 
general. Interpolation is awkward. A rather large computer is 
required to operate this system efficiently. 

c) Horphologic lines and points 

Significant terrain points (e.g., peaks, troughs) and 
terrain break lines (e.g., ridges, rivers, road edges) are measu- 
red, which gives a more or less random point pattern. 

Advantages : A minimum number of points are measured 

Valuable in delineating man made features. 

Disadvantages : The point selection would have to be carried 

out by a trained engineer or photogramiiiotr ist . 
The sampling cannot be automated. 

Interpolation is elaborate. 

d) Heterogeneous grids 

B.Hakarvic developed a method (progressive sampling) 
where a square is used as pattern of registration but the height 
is measured only at those intersection points where it is 
necessary .This very promising method is described in the next 


chapter . 
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3.3 Sampling denaity 

It is possible to achieve any accuracy of the desired output 
of a DTM systein(i.e. the reconstructed terrain) by varying the 
density of sampling, no matter what pattern of registration is 
used. Any increase in accuracy will .however , mean that the cost 
of acquisition will also increase. Normally for flat terrain the 

V 

data points are collected at a distance, for undulated the data 
points are collected closely. 

3.4 Data processing 

The problem of deriving from the stored terrain data 
( e.g. contours or heights in a regular grid, or profiles) some 
special information (e.g. data for hill shading, volumeti-i c 
data, contours .... ) is first of all one of interpolation. 

A number of different methods are known and used in practice 
to determine the height for any intermediate point from the given 
X,Y,Z coordinates. They range from linear interpolation to quite 
sophist icated concepts like linear least squares 

interpolation. Detailed description of various interpolation 
methods are described in next chapter. 
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CHAPTER IV 


MATHEMATICAL FORMULATION AND ALGORITHMS 


4.1 Introduction 


In this ch«ipt*r thtt description of various slftorithns of 
aathenaticnl procedures that are used in this is work are 
presented. They are 

(1) weiflhted averaae method 

(2) Interpolation methods 

(3) Progressive sampling method 

The first step in DTM is to convert discrete data into 
gridded pattern. Size of a grid is usually taken as 20 mts. 
Normally in DTM photographs or maps of 1:10000 scale are used 
for data collection. So the size of a grid on a photograph or on 
a map will be 2mm which can be easily measured. If we take 10 mts 
grid more points are generatr-d and storage becomes difficult. 

4.2 Weighted average method Algorithm 

(1) Gridding is the process of determining the values of the 
surface at a set of locations that are arranged in a regular 
pattern, usually a square from irregularly located control points 
where the values of the surface are known. 

(2) As the first step control points must be sorted according to ' 
their geographic coordinates. Normally four control points are 
used to estimate elevation at grid node as four points are 
assumed to be in four directions. 

(3) Fig 4.1Ca) shows a series of observations on map. Each point 
is characterized by its X coordinate, Y coordinate and Z 




(a) Oriof,aa.l set of irregularly spaced points on a. wap. 

(b) Regularly spaced grid netta>ork, with nodes whose values 
are to be calculated* 

(c) Location of four nearest control poiimts to one gyld 

node* 

Cd) Cowplftt® grid with estSniated elevations at every node. 
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coordinate or elevation. The observations may be identified by 

numbering them sequentially as they are read in from 1 to 1. 

Therefore an original data point I has coordinates Xi,Yi,and Zi. 

(4) In Fig. 4. 1b, we have superimposed a regular grid of nodes on 

the map. These grid nodes are numbered sequentially, from 1 to 

k. Grid node k has coordinates Xj^, Y|( and has an estimated 

elevation of Zj^. 

(5) Now the next step is to search four nearest data points for 
each grid point as shown in Fig. 4. 1c. 

(A) After searching four nearest data points to grid node k the 
distance D^i^. from observation point i to grid node k is found by 
Pythogorous equation: 

Dik = C(Xk - Xi)2 + (Yk - 

Having found the distances Dik to four nearest data points 
we now can estimate the grid node elevation Zk from these. The 
estimate is 

Zk = <Zi/Dik) / (1/Dik> 

(7) In a similar manner we can evaluate the remaining grid nodes 
on the map. The computed grid with all values of Zk posted is 
shown in Fig. 4. Id. 

4.3 Interp olat i o n Procedure 

Using the method of weighted average method the elevations at 
four corners of each grid can be calculated. In progressive 
sampling to locate new point we need a grid with nine points 
called "Zero run" and from X,Y coordinates of the new point we 
should calculate Z value. To solve this problem the best way is 
to express Z as a function of X,Y so that structure required for 
Zero run as well as for further interpolation Z values can be 
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Fig. 4.2 GRAPHICAL DEPICTION OF 
LINEAR INTERPOLATION 
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There are a number of different interpolation procedures to 
determine the height for any intermediate point from the given 
XyY,Z coordinates. They are 
(I) Linear Interpolation 
<2> Least Square Interpolation 

4.3.1(a) Linear Interpolation 

The simplest form of interpolation is to connect two data 
points with a straight line. The technique is called linear 
interpolation and is depicted graphically in Fig. 4. 2. 

Using similar triangles 

Cfi(X) - f (Xo)3/(X-Xo) = Cf<X^) - f (Xj,)3/(X-|-Xo) 
which can be rearranged to yield 

f^CX) = f<Xo) + fCf(Xi) - f (Xo>3/(Xi-Xo>> (X - Xq) 

In general, the smaller the interval between the data 
points, the better the approximation. 

4.3.1(b) Bilinear Interpolation 

Here it is assumed that Z varies linearly in both X and Y 
directions. 

A bilinear polynomial Z *= A + BX + CY + DXY would 
suit a undulated terrain. This is illustrated in Fig. 4. 3. 

The bilinear constants are 

A = Zi B = (Z2-Zi)/deltaX 

C » (Z3-Zi)/deltaY 

and D = (Zi+Z 4 -Z 2 -Z 3 )/(deltaX*deltaY) 

Bilinear constants are calculated and they are used for 
further interpolation in progressive sampling and Earthwork 


calculation 
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4.3.g. Least square Interpolation 

This method was developed by Legendre over 150 years ago and 
this has been used in Statistics for many years. Earth, assumed 
to be a plane surface was represented by mathematical model 

Z = f(X,Y) = L + MX + NY 

L,M,N are the parameters and (X,Y,Z} represent any point on 
the plane surface. 

Graphical depiction of multiple linear regression where Z is 
a linear function of X and Y is shown in Fig. 4. 4. 

The best values of coefficients are determined by setting up 
the sum of the squares of residuals 

Sr = r <Zi - L - MXi - NYi)^ 
and differentiating with respect to each coefficient 
SSf./ 6L = -2 Z (Zi - L - MXi “ NYi> “ 

tfSp/dM « -2 E Xi(Zi - L - MXi “* ^Yi) (4.1) 

dSr/«N = -2 E Yi(Zi - L - MXi “ NYi) 

The coefficients yielding the minimum sum of the squares of 
the residuals are obtained by setting the partial derivatives 
equal to zero and expressing Eq.(4.1) as a set of simultaneous 
linear equations. 

nL + M E Xi + N E Yi = E Zi 
LE Xi + M E xf + N E XiYi « E XiZi 
L E Yi + M E XiYi + N E Yi « E YiZi 

For flat terrain where height variation is less than 10%, 
this model can be adopted. These constants are used for further 
interpolation. 
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Fig.4.4 GRAPHICAL DEPICTION OF 
MULTIPLE' LINEAR REGRESSION 
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4.3.2. (b) Collocation 

The method of collocation, sometimes also known as advanced 
least squares interpolation has been extensively applied to 
solution problems related to geodesy and photogrammetry. The 
back bone of the method is so called covariance function, which 
takes into account the effect of correlation between 
observations. 

The basic assumption made here is that for a pair of points 
on the surface of earth, the pair of heights are correlated 
according to the distance between points. This means that the 
distance and direction are involved in the function which is 
ultimately to be used for derivation of height differences. 
Thus, if we have a set of given points with heights 
h1,h2,h3, . . . . .hn, then height *h* of an unknown point can be 
written in function form as 

h* f(r1,r2, rn, ©1,02, ©n) 

where r1,r2.....rn and ©1 ,62, ...... .On are the distances and 

directions of the given points from unknown points (height only). 

This method is elaborately described in references quoted 
under C14,16,173 


* 
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4,3.3 Statistical Analysis 


Statistical analysis is carried out to express the accuracy 
of the model. 

As a first step variance covariance matrix for regression 
coefficients is computed. 

Variance is computed by dividing corrected sum of squares by 
degree of freedom. Variance is defined as the average square 
deviation from mean. 

covariance is the Joint variation of two variable about 
their common mean. This measures the distribution of values 
around common mean. 

Covariance is calculated by dividing corrected sum of 
products by degree of freedom. 

Secondly correlation matrix for regression coefficients is 
computed. To estimate the degree of interrelation between 
variables correlation coefficient *r* is calculated. Correlation 
is the ratio of covariance of two variables to the product of 
their standard deviations. 

Posterial' variance ((r*^o > is calculated to express absolute 
accuracy. 

1/E 

O" *o » C (Sum of squares of residuals ) /(n—l )□ 

4.4 Progressive Sampling 

Progressive sampling of discrete points arranged in a 
regular grid of varying density provides a powerful tool for the 
acquisition of digital terrain data. The sampling process is 
started by measuring a low density point grid. The corresponding 
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data are analyzed to locate the positions of new points to be 
sampled in the next run. 

The algorithm and logic of the method are explained below. 
Algorithm 

(1) For organizational reasons and to reduce the required 
computer storage the entire area is divided into grids. 

(2) The heights at four corners are known to us from the output 
of weighted average method. The heights at nine points are 
calculated. 

This structure is called zeroth run and is shown in fig4.5. 

(3) Progressive sampling is supposed to be carried successively 
for all grids in an ordered manner. For transition from one grid 
to next, the boundary data should be kept temporarily in the 
computer storage. 

(4) In numerical processing the first step of the analysis is 
computation of height differences between the adjacent points 
along each row and column 

Zi.j-8 = 

Zi-8.j = 

(IS) The next step is to calculate second differences separately 
for the rows and columns. This can be formulated algebraically 
for the point in the following way 

Z-I .25 (k.i-D - 3^1.9 (k.i) = 6 (r) for a row 

^2S.1 1.i) — Z 9 ,-| (k.i) = 6 Z-|,.| (c) for a column 

The heights sampled in the following runs are processed 
similarly. When disregarding random terrain fluctuations. the 
second differences carry information about terrain curvature. 
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When curvature is substantial, it is desirable te 
sampling density. 


increase 


the 


(6) For the selection of new points the following criteria may 
be applied. If the 5 Zi.j.LE.tf th no further sampling is needed 
If the second difference th exceeds a certain threshold the 
sampling density should be increased. The values of threshold 
for various types of terrains are taken from reference (4 ). 

<7> After the central points to be sampled are determined, the 
surrounding peripheral points should be considered for the next 
run. Steps (1) to (6) are to be repeated for all runs. 

<S) The low density grid with 3x3 points is considered first 
and IS densified upto a maximum of 33x33 in successive runs. 

(9) The same procedure is repeated for all other grids. 

(10) The flow chart of computer operations involved is given in 

next page. 
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CHAPTER S 

VOLUMETRIC COMPUTATION 

5.1 INTRODUCTION 

In this chapter the application of DTM to volumetric 
computation is emphasized. The methods mentioned in chapter two 
are elaborated here. 

■ The ultimate purpose of field surveys is to obtain certain 
numerical values to represent quantities such as areas and 
volumes. Surveyors are often called upon to measure volumes for 
engineering projects viz.. highway design, water resources 
structures construction and building construction. Volumetric 
determination is an important aspect in practical surveying. 
Earth work operations involve the determination of the volumes 
of materials which must be excavated or embanked on engineering 
projects to bring the ground surface to predetermined grade. The 
field work involves the measurements of the dimensions of various 
geometrical solids which make up volumes, the setting of grade 
stakes, and the keeping of field notes. The office work involves 
the computations of the measured volumes and determination of 
most economical measures of performing work. Because of the 
inherent accuracy of modern topographic maps of large scale 
produced by photogrammetr ic methods much of the field work is' 


reduced . 
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Volumetric determination by point grids of digital terrain 
model data offers the advantage of great flexibility in route 
location. The main reason for this is it has become possible to 
collect the point data from photographs in automatic modes of 
equal spacing using analytical and semi automatic 
photogrammetric instruments. 

5.2 Volumetric computation using DTM : 

The basic task of a highway engineer is to select a proper 
alignment between given two termini. For that he has to study 
(i) road alignment search procedure (ii> sight distance 
computation and (iii) earthwork computation. Then optimal route 
(s) will be decided. Surveyor plays an important role in earth 
work calculations. 

Volumetric determinations are essential in highway 
design phase. More commonly, actual ground measurements are 
obtained for such determinations. Application of modern 
techniques of photogrammet ry to highway design in developing 
countries is still more or less in an initial stage due to lack 
of knowhow. Research on this field is therefore limited. Use of 
photographs provides us additional information about types of 
soils and geological characteristics of terrain under 
observation. 
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The terrain surface is essentially defined by complete 
numerical study in X,Y, and Z coordinates of very large number 
of points distributed over the strip selected for the projects. 

In developed countries, numerical photogrammet ry and 

automatic computations are being more and more adopted for 
alignment and designing motorways. The X,Y,Z values taken as 
input to computer with a suitable mathematical processing model 
to obtain the volumes of cut and fill. 

For the present study cross section area method is used for 
volume computation. Different cross sections are shown in fig 
5.1. The cross section shown fig 5.1 (a) is adopted for volume 
calculation for flat terrain. For undulated terrain the cross 
sections shown in fig 5. 1 ( b) ,5. 1 ( c ) can be used for volume 
calculation. 

The logic of the computer program is shown in the 
flow chart. The results are presented in chapter 6. 




Fie. 5.1. 
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START 3 

r 4—=^ 

READ DTM, SIDE SLOPES, WIDTH OF ROAD 
TOTAL POINTS, CENTER LINE COORDINATES 
LENGTH OF CROSS SECTION. 


^ 

DO 100 1= 1, TPS 


1 COMPUTE ZDTM, ZCOR, AREA(I) 





|100 CONTINUE 



— . i. 

DO 110 K= 1, TPS-1 


: ::: 4 . . 

VOL(K)« 0.5*<AREA(K>+ AREA{K+1 ) - 


110 CONTINUE 



STOP 


FLOW CHART FOR VOLUME COMPUTATION 
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CHAPTER VI 
RESULTS 

6*1 Analysis 

In this chapter results are presented. Fig 6.1. shows the 
locations of data points collected for flat terrain. The K,y and 
2 values of the data points are given in Table 6.1.1. The 
locations of points to form the model are shown in fig 6.2. 
Totally twenty points are taken to solve the model Z= L+ MX + NY, 
so that there exists a degree of freedom of 17.'irhe constants L, 
M and N are computed using STAT package. PosteriaV variance is 
computed to indicate absolute accuracy. The results are shown in 
Table 6.1.2. 


Fig 6.3 shows the locations of data points collected for 
undulated terrain. The X, Y and Z values of the data points are 
given in Table 6.1.3. To analyse undulated terrain entire area is 
divided into five parts, as shown in fig 6.4. Five sets of each. 
Variance- Covariance matrices, correlation matrices for 
regression coefficients and posterial' variances are computed. 
They are shown in Table, 6. 1.4 to Table 6.1.8. 



54 


6.2 Siaulated terrain generation: 

An oblique view of a surface provides a visually 
effective representation of slopes and peaks. Displaying a space 
curve or surface is not an easy task. The simplest method of 
displaying curves and surfaces is wire-frame technique. To 

display a curve, its coordinates are evaluated for many closely 
spaced points and these points are connected by short straight 
line segments. A complete surface is displayed by generating a 
mesh of such curves holding one parameter constant at a time. 

Geometric transformations play an important role in 

generating images of three dimensional scenes. The important 

geometric transformations are translation, rotation and scaling. 
Then perspective transformations are carriedout to display the 
three dimensional surface. 

In the present investigation three dimensional terrain 
is simulated using PLOTIT software on IBM PC-AT with color 
monitor. The sequence of computer operations involved are shown 


in flow chart 



START 


ENTER PLOTIT 


INPUT DISCRETE X,Y,Z VALUES, NO. OF COLUMS 
AND NO. OF ROWS. 


INPUT DEBIGW PARAMETERS, POLYNOMIAL GOVERNING 
THE SURFACE, VIEW ANGLE, BASE HEIGHT, 
PERSPECTIVE FACTOR. 


INTERPOLATION 

TRANSFORMATIONS 

DISPLAY THE SURFACE 

STOP 

FLOW CHART FOR COMPUTER OPERATIONS IN 
SURFACE MODELING. 
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6.3 Volumetric Computation 

Volumes are computed fay using end area method. Three 
different routes between two given control points passing through 
different obligatery points are considered for volume 
calculation. 

The perameters for volume calculation are digital terrain 
midel, side slopes, width of the road and length of cross 
section. , 

For the present study flat terrain is considered. The other 
design parameters are 

Side slopes =1!1.5 

Width of road =12 mts 

Length of cross section = 40 mts 

Tables 6.1.9 to 6.1.11 shows the results. Route I is the 
optimal route between two given control points with minimum earth 


work 
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6.3 Mass diagram: 

Mass diagram or mass curve is a curve plotted on a 
distance base, the ordinate at any point of which represents the 
algebraic sum up to that point of the volumes of cuttings and 
fills from start of the earthwork or from any orbitrary point. In 
obtaining algrbraic sum the convention is to consider fills 
negative and cuts positive. The positive total volumes are 
plotted above and the negative total volumes are plotted below 
the base line. 

Uses of mass curve: 

With the sign convention adopted upward slope of the 
cuve in the direction of algrbraic summation indicates excavation 
as downward slope indicates embankment. 

Between any two points in which the curve cuts the base 
line the volume of excavation equals that of embankment, since 
the algebraic sum of the quantities between such points is zero. 

The vertical distance between a maximum point and the 
next forward minimum point represents the whole volume of an 
embankment, that between a minimum and next forward maximum point 
the whole volume is cutting. 
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la Analysis of flat terrain: 

The Variance- Covariance matrix for 

Regression Coefficients 

Coeff . 

L 

M 

N 

L 

0.070 

-1 .30E-04 

-1.77E-04 

M 

-1.30E-04 

6.50E-07 

o 

• 

o 

N 

-1.77E-04 

o 

• 

o 

1.15E-06 

Correlation 

matrix for 

Regression Coefficients* 

Coeff. 

L 

M 

N 

L 

1.00 

-0.607 

-0.620 

M 

-0.607 

1.00 

o 

m 

o 

N 

-6.20 

0.0 

1.0 


Regression Coef f i c i ents s 

L = 204.551 
M = -8.33E-05 


N = -0.001 



Table 6.1.2: MODEL Z =L+MX+NY 


Z observed 

Z fitted 

204.50 

204.53 

205.30 

204.49 

204.20 

204.51 

205.40 

204.22 

204.10 

204.37 

205.30 

204.45 

205.30 

204.24 

205.30 

204.23 

204.40 

204.42 

204.00 

204.35 


PosteriaV variance = 0.7308729 
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Table 6.1.2(a) 

MODEL Z=A+BX+CY+DXY 

Z observed Z fitted 

207.2070000 
209.2020000 
209.9230000 
208.9430000 
209.5870000 
209.5170000 
207,4410000 
209.6150000 
209.9230000 
209.1530000 

Posteriajr Variance = 4.3203000 


204.5000000 

205.3000000 

204.2000000 

205.4000000 
204.1000000 
205.3000000 
205.3000000 
205.3000000 

204.4000000 
204.0000000 







INPUT DAT4 FOR UnDULATFO IERRATIN 
table 6.1.3, 


il-5 

51 

198 5 

bBa^a- 

48,0 

1989 

157 

22.0 

1994 

?22..> 

58,0 

1990 

2 8 B , r" 

47.5 

2002 

3p5,fc 

4 8.0 

2 000 

7 0.-. 

105.0 

2005 


88.0 

1987 

3 47,5 

92.0 

2002 

3i>, » 

138, 

2000 

192 

129 

1996 

7l;U 

168,0 

1 990 


159.0 

1982 

363. .V 

168,0 

1987.5 

185. 

198,0 

1994 

■2 60 

208.5 

1990 


227.5 

1997,3 

, 

S20.r^ 

228.0 

199 0 

no,& 

250.0 

1990 

380 

231 

2001 

362 

260 

2003 


306 

2 000 

121 

312 

1986.4 

'220 

301 

2003.5 

bi 

384 

2002 

^53 

312 

1999 

3 5.':? 

383 

1987.2 

bl, a 

365 

1989 

110 

391 

1987 

16 4..T. 

390 

1989 

'2I3.U,'' 

385 

1994.3 

332 

383 

2001 

362 

386 

1997 



Fifi 6.4 0 LOCATIONS OF DATA POINTS TO FORM MODEL 

Ei LOCATIONS OF POINTS TO TEST MODEL 
SCALE 1 cni = 20 mts. 
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6.1.1b Analysis of undulated terrain: 

SET I: 

The Variance- Covariance matrix for Regression coefficients 


Coef f . 

A 

B 

C 

D 

A 

1.111 

-0.056 

-0.033 

0.002 

B 

-0.056 

0.006 

0.002 

-1.67E 

C 

-0.033 

0.002 

0.002 

-8.33E 

D 

0.002 

-1.67E-04 ■ 

-8.33E- 

-05 8.33E 

Correlation 

matrix for 

Regression 

Coefficients! 

Coeff . 

A 

B 

C 

D 

A 

1.000 

-0.707 

-0.775 

0.548 

B 

-0.707 

1.000 

0.548 

-0.775 

C 

-0.775 

0.548 

1.000 

-0.707 

D 

0.548 

-0.775 

-0.707 

1.000 

Regression 

Coef f i cients s 




A 

= 1999.333 




B 

= 0.183 




C 

= 0.100 




D = 2.23E-09 
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SET II! 

The Variance- Covariance matrix for Regression coefficients! 


Coeff . 

A 

B 

- - - - — ' 

C 

D 

A 

6.267 

-0.313 

-0. 188 

0.009 

B 

-0.313 

0.031 

0.009 

-9.40E-04 

C 

-0. 188 

0.009 

0.009 

-4.70E-04 

D 

0.009 

-9.40E-04 - 

-4.70E- 

-04 4.70E-05 

Correlation 

matrix for 

Regression 

Coefficients! 

Coeff. 

A 

B 

C 

D 

A 

1.000 

-0.707 

-0.775 

0.548 

B 

-0.707 

1.000 

0.548 

-0.775 

C 

-0.775 

0.548 

1.000 

-0.707 

D 

0.548 

-0.775 

-0.707 

1.000 

Regression 

Coef f i cients 

A 

B 

C 

• 

> • 

= 1990.000 

= -0.204 

= 2.38E-08 




D = 0.006 
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SET III: 

The Variance- Covariance matriK for Regression coefficients: 


Coeff . 

A 

B 

C 

D 

A 

20.851 

-1.043 

-0.626 

0.031 

B 

-1.043 

0.104 

0.031 

-0.003 

C 

-0.626 

0.031 

0.031 

-0.002 

D 

0.031 

-0.003 

-0.002 

1.56E 

Correlation 

matrix for 

Regression 

Coefficients s 

Coeff. 

A 

B 

C 

D 

A 

1.000 

-0.707 

-0.775 

0.548 

B 

-0.707 

1.000 

0.548 

-0.775 

C 

-0.775 

0.548 

1.000 

-0.707 

D 

0.548 

-0.775 

-0.707 

1.000 


Regression Coefficients: 

A = 1992.000 
B = -0.396 
C = 1.19E-07 


D = 0.027 



SET IVs 

The Variance- Covariance matrix for Regression coefficients! 


Coeff . 

A 

B 

C 

D 

A 

22.778 

-1.139 

-0.683 

0.034 

B 

-1 . 139 

0.114 

0.034 

-0.003 

C 

-0.683 

0.034 

0.034 

-0.002 

D 

0.034 

-0.003 

-0.002 

1.7 IE-04 

Correlation 

matrix for 

Regress! 

on Coefficients! 

Coef f . 

A 

B 

C 

D 

A 

1.000 

-0.707 

-0.775 

0.548 

B 

-0.707 

1.000 

0.548 

-0.775 

C 

-0.775 

0.548 

1.000 

-0.707 

D 

0.548 

-0.775 

-0.707 

1.000 


Regression Coefficients! 

A = 2000.333 
B = 0.233 
C = 0.12S 
D = -0.025 
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SET V: 

The Variance- Covariance matrix for Regression coefficients! 


Coeff. 

A 

B 

C 

D 

A 

9.701 

-0.485 

-0.291 

0.015 

B 

-0.485 

0.049 

0.015 

-0.001 

C 

-0.291 

0.015 

0.015 

-7.2SE-04 

D 

0.015 

-0.001 

-7.2aE-04 

7.2aE-05 


Correlation matrix for Regression Coefficients! 
Coeff. A B C D 


A 

1.000 

-0.707 -0.775 

0.548 

B 

-0.707 

1.000 0.548 

-0.775 

C 

-0.775 

0.548 1.000 

-0.707 

D 

0.548 

-0.775 -0.707 

1.000 

Regression 

Coeff ic 

ients! 

A = 2004.450 

B = -0.039 

C = -0.217 



D = 0.008 



Table 6.1.4 


Z = A+BX+CY+DXY 


Z observed Z fitted 


1998.00 

1999.33 

2002.80 

200.1.33 

2012.00 

2013.33 

2003.40 

2003.33 

2004.80 

2005.33 

2008.50 

2007.33 

2007.00 

2007.43 

2014.00 

2014.90 

2010.50 

2011.33 

2012.80 

2013.33 


PosteriaJ' variance 


1.0657390 
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Table 6.1.4(a) 


Model Z = L+MX+NY 


Z observed 


Z fitted 


1998.0000000 

2000.3330000 

2002.8000000 

2004.4330000 

2012.0000000 

2010.0300000 

2003.4000000 

2000.4000000 

2004.8000000 

2006.7960000 

2008.5000000 

2010.6540000 

2007.0000000 

2005.4320000 

2014.0000000 

2016.3330000 

2010.5000000 

2012.4000000 

2012.8000000 

2010.7980000 




PosteriaV Varaince 


2.23 




Table 6.1.5 


HODEL Z = A+BX+CY+DXY 


Z observed 


Z fitted 


1973.40 

1973.68 

1983.90 

1983.28 

1985.50 

1986.64 

1996.00 

1994.32 

1992.10 

1993.12 

2010.80 

2012.08 

2006.00 

2006.56 

2018.00 

2017.60 

1999.00 

1997.92 

2020.00 

2021.68 


Posteriay variance =* 1.1414550 
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Table 6.1.5(a) 


Model Z = L+MX+NY 


Z observed 

1973.4000000 

1983.9000000 

1985.5000000 

1996.0000000 
1992.1000000 
2010.8000000 

2006.0000000 

2018.0000000 

1999.0000000 

2020.0000000 


Z fitted 

1975.6440000 

1985.5000000 

1985.4990000 

1993.9200000 

1990.1330000 
2012.8330000 
2004.5350000 
2016.4000000 

2000.1330000 
2018.4950000 


Posterial' Varaince = 1.96 



Table 6 . 1 • 6 


MODEL Z = A+BX+CY+DXY 


Z observed 

Z fitted 

1975.90 

1976.16 

2010.00 

2009.28 

2168.40 

2168.40 

2015.90 

2015.04 

2015.90 

2016.48 

1969.00 

1968.24 

1936.80 

1936.56 

2239.00 

2238.96 

2153.90 

2154.90 

2116. 10 

2115.12 


Posteria^ variance 


0.6879366 



Table 6.1.6. (a) 


Model Z = L+MX+NY 


Z observed 

Z fitted 

1975.9000000 

1973.2000000 

2010.0000000 

2012.5340000 

2168.4000000 

2170.5330000 

2015.9000000 

2017.8950000 

2015.9000000 

2016.9350000 

1969.0000000 

1971.4950000 

1936.8000000 

1938.5330000 

2239.0000000 

2241.5350000 

2153.9000000 

2156.7950000 

2116.1000000 

2119.2940000 


Poslerial' Varaince = 2.48 





Table 6.1.7 

MODEL Z = A+BX+Cy+DXY 


Z observed 

Z fitted 

1990.90 

1992.15 

1968.00 

1967.49 

1940.00 

1939.65 

1981.50 

1981.47 

2028.00 

2028.29 

1929.30 

1928.63 

1955.00 

1954.77 

1881.80 

1882.61 

1761.90 

1762.95 

1834.00 

1833.63 


Posterialf variance 


0.7013458 






7 9 • 


Table 6.1.7(a) 


Model Z = L+MX+NY 


Z observed 

Z fitted 

1990.9000000 

1994.8950000 

1968.0000000 

1964.5350000 

1940.0000000 

1936.6950000 

1981.5000000 

1983.4950000 

2028.0000000 

2030.5330000 

1929.3000000 

1923.5650000 

1955.0000000 

1960.4550000 

1881.8000000 

1879.7990000 

1761.9000000 

1760.4550000 

1834.0000000 

1835.2350000 


PosteriaV Varaince 


3.65 
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Table 6.1.8 


MODEL Z =A+BX+CY+DXY 


Z observed 

Z fitted 

2001.00 

2000.25 

2005.00 

2004.95 

2003 . 00 

2002.11 

2011.90 

2012.21 

2019.00 

2018.25 

2018.00 

2017.89 

2052.50 

2053.17 

2020.90 

2021.89 

1997.00 

1996.83 

2130.00 

2127.34 


PosteriaJf Variance 


1.0832610 
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Table 6.1.8 (a) 


Model Z = L+MX+NY 


Z observed 

Z fitted 

2001.0000000 

2000.3450000 

2005.0000000 

2003.4550000 

2003.0000000 

2001.8980000 

2111.9000000 

2010.9000000 

2019.0000000 

2020.7430000 

2018.0000000 

2023.9980000 

2052.5000000 

2055.5430000 

2020.9000000 

2022.9500000 

1997.0000000 

1998.3330000 

2130.0000000 

2132.4550000 


Posteriay Varaince = 2.5 
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ROUTE I! 
TABLE 6.1.9: 


SECTION NO. 

PILL (cubic 

CUT (cubic at«) 

1 


7.178 

2. 


5.587 

3. 


3.991 

4. 


2.393 

5. 


0.798 

6. 

0.803 


7. 

2.404 


8 . 

4.003 


9. 

5.607 


10. 

7.213 


TOTAL VOLUHE « 39.174 cubic site. 
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ROUTE II: 
TABLE 6.1.10: 


SECTION NO. 

FILL (cubic otta) 

CUT (cubic nt«) 

1. 


12.600 

2. 


6.234 

3. 


2.243 

4. 


0.644 

5. 

3.356 


6. 

7.364 


7. 

6.56 


8 . 

5.758 


9. 

4.954 


10. 

1.748 



TOTAL VOLUME « 51.36 (cubic mts) 
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ROUTE III: 
TABLE 6.1.11. 


SECTION NO. 

FILL (cubic nta) 

CUT(cubic mts). 

1. 


5.871 

2. 


9.067 

3. 


9.861 

4. 


10.656 

B. 


6.676 

6. 


0.286 

7. 

3.717 


8. 

7.726 


9. 

11.743 


10. 

13.354 



TOTAL VOLUME = 78.914 (cubic nts) 
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CHAPTER VII 

CONCLUSIONS AND SUGGESTIONS FOR FURTHER RESEARCH 

In the work presented an attempt has been made to 
explain the methods to generate the digital terrain model and 
using it for volumetric computation for highway alignement with 
the use of digital computers. The following are concluded from 
the present study. 

7.1 Conclusions! 

We can generate DTM for truely plane areas with 
sufficient absolute accuracy(0 o= 0.73m) with the model 
Z= L + M X + NY. 

Similar accuracies may be possible for larger areas and 
cannot be commented upon since not tested. 

WE can generate DTM for undulated terrain by sudividing 
the area into five parts with the model 

Z=' A ♦ B X-! C Y + D X Y. ' 

Progressive sampling provides a powerful tool to 
generate a cartographic database. 

Computer grapfiics is a key element for decig'n And as€ of 
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geographical information system. Simulated three dimensional view 
of terrain "fiows that graphical product and numerical output are 
in agreement with each other. 

Volumetric computation through DTM provides a powerful 
tool for alternate alignment selection. 

The rapid development in the ability to handle terrain 
data in a complete digital form holds forth the promise of 
reducing the drudgery of cartographic operations of providing a 
wide variety of data interactions, and of reducing time and cost 
so that managers and decision makers will know how to make 
maximum utility of resources of world. 

DTM provides a base for geographical information system. 

7.2 Suggest f o 1 1 s f or future works 


In the present study analysis of flat and undulated 
terrain are carried out. A suitable mathematical model can be 
gerenated for hilly terrain. 

For volumetric computation in undulated and hilly areas 
the length of cross section should be taken at small intervals. 
Additional corrections such as curvature corrections should be 
considered. 

Spot heights from progressive sampling data can be 
effectively used for corridor study. 

The DTM can be used to rectify satellite born imagery 
as explained in ref erenceCHT . 

Contouring can be done by using DTM. 
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APPENDIX 



t3«o0i«tric Trans fornat ions 
Translation 

The transformation which translates a point (x,y,z) to a new 
point (x’ ,y’z’ ) is 

[X’ y’z' 1] = [X y z 1] 1 0 0 0 

0 10 0 
0 0 10 
Tx Ty Tz 1 

where Tx,Ty, and Tz are components of translation in x,y and z 
directions respectively. 

Rotation 

In three dimensions, it is helpful to device transformations 
for rotation about each of the three coordinate axes. Rotation 
about the Z coordinate axis through angle theta is 

[x’ y’ z’ 1] = [x’ y’ z’ 1] cosO -sinO 0 0 

sinB cosO 0 0 

0 0 10 

0 0 0 1 

The rotation angle 6 is measured clockwise about the origin 
when looking at the origin from a point on the +Z axis. Similarly 
rotation about X and Y axes cna be written. 



Scaling 


A acaling tranaformatl 
in each coordinate directions 

tx’ y» z' 1] = [X y z i] 


on can be used 

separately 

Sx 0 0 

0 Sy 0 

0 0 Sz 

0 0 0 


to scale 


0 

0 

0 

1 


dimensions 


Perspective Transformation 

A perspective display can be generated by simply 
projecting each point of an object onto the plane of display 
screen as shown fig. The coordinates (Xs, Ys) of the projected 
image of point P measured in eye coordinates CXe,Ye,Ze) are 

Xs = (D*Xe/S*2e) 

Ya = (D*Ye/S*2e) 



>'E THREE PRIMITIVE THREE DIMENSIONAL ROTATIONS 



P(Xc ZO 




ti 


C 


C 


C 

80 


C 


K02 

C 


104 


103 

C 

10 


101 

100 

81 

U 


PRtlGF.AMMF T3 CONVERT A OISOkETF OA^A IM GRTOOeO PATTERN 

Py O.V.*lAGARAJtJ 

*^*4********-r**^'^-.^*^*:^*itf.t^f.t’^^,f:i^*^^*it$*$^^*t* ********** 

DEOL.AHATIOM 
real A(100,3’> 
real Bf250,2501 

real MIN)C,MINlC,MAXX,MAXy - OFLTAX.DELTAy 

real X, y-SMAT.L 

REAL DlSrciOO) 

integer «AXPS,MAXRS,maXCB 

OPEN 

0FEN(UMIT=71 .FTLEs'UPIN.nfiT") 

OPEiaUMIT = ?4,FXLF='UKOOT,UAT') 
npRLf UMITz?3,FILE='El4E,DAT') 
read STATEMEfaXS 
READf2},») MAXPS, MAXES, MAXES 
READf21,T) MINX, Miwy^HAXX.MAXy , small 
READ( 21 ,») C(A(t,J), J=1,35,I=1,MAXPS) 
DEf-rAX=(HAXX-MlHX)/CMAXCS-l ) 
DELrAy=CMAXy-HTNY)/(MAXRS-l) 

MAIN LOUP 
WRITE (74,80) 

FORMATCieX, '’X*,9X, 'y',7X,'Z'///) 


y - MiNY 

DO 100 1 = 1,MAXRS 
XsMiPX 

r>0 H)1 3 St, M Axes 

calculate DIST’t't^ BFTWEEN CPRRENT GRIDPDiNi AND ALL DATA POINTS 
Do 107 K =i,MAXPS 

n s fX-AC K,1 ) )4t2+(Y-A( K,2))*47 
D1SX(K)=:S0R1 f D1 
CDNTTHOE 

find the nearest 4nATA pniNTS and calulate spms 

S1=0.0 
52—0 .0 

DU 103 K=l,4 
ir=j 

DO 104 L=2, MAXES 

TFfDISTtLl.LT.DiSTflC)) IC=L 
CUNTTNHE 

IF CDlST(ir).LT, SMALT.) GD To 10 
.S1=S!*ACIC,3)/D1ST(IC) 
s2=S7fl .o/nisirir) 

DlSTflD = 9790999,90 

continue 

calculate grid point and store in matrix 


Z-S1/S2 
Z =A(TC,3) 

WRTTF(73,t') 7. 

TF(J.EO.l) WRITE(2^,B1) T,X,Y,Z 
XFCj.NF.n WRiTE(24,P2) X,y;z 
TF fJ.En.STI WRITE(74,83) 

X=X+DELTAX 

CONTINUF 

YsY+DELTAY 

continue 

FORMAT(10X,I2,2X,FR,3,3X,F6 1,3X,F6.'1) 
FORMAT (14X,FP. 3, 1X|F6,1,5x,F 6,1) 
FORMATf 13X,2RC'-')) 

STOP 

end 



lt?U 03 f 

tsoOiV 

DOb’' 
ODD? 
ttbOBf" 
»l>09i' 
PUl tiT 
t?ci lt.■ 

a&i 3^ 
II& 14 ' 
5' 

tJUl b v- 
tiul 7 t 
BOJ 8 
1*{>J 9‘ 
tlD 2 y:' 
002 i‘ 
002 2 
0023* 
0024 
W025*'' 
O 02 br 

002 7' 
0028=' • 
00 2 9? 
0 03 v:“. 
Oo3l> 
tl0 3 2< 
0033^ 
00 3 4-: 
Oo3b''- 
oa3b'- 

003 7 
I? 0 3 8 i' 


s# 

c 

c 

c 

c 


c 

77 

88 

99 


30 

50 


tntege;i^'dp5’s,l,m,n 

REAL SCRR.FUM.DEN 

norl!! =2I,FILF = ' T fjSl A r . D A'j- ' ) 

reaS ""‘'"•“r') 

dES‘==[)pi 5 !i '=''"tO,i=i,opfsi 

MAIN PRDGRAM»-«f; 

Du 77 isj^nprs 

D3 88 1 “*= l.npTs 
PCT) =.(Zn(T)-ZF(I)U»2 
fQDl'TNnE 
PS =f--.i 

DD 99 I = I,DPTS 
SS =,5S +S(T) 
rOr.’l'lNf’E 
FUR sSP/DFf! 
vPERk =P3RTfFtlN) 

'^'RIIF C24,30) ( ZO C 1 ) - T = 1 - DpTS ) 

WRII'F (24;3u) ahl) IsI'dptS) 
FURHM' (5j(, 3FJ0.8) ' i.RPTSJ 

WRTI'F (24,50) KERR 
F3RMAJ CfSX, FH'.2) 

Sinp 

FND 



c 

f 

c 

c 

c 

c 


c 


c 

c 


c 

c 

c 


kU 

c 




c 


c 


bi> 

c 

7 0 


m*'i^***t^*i'***^**^*^***^^*^*^****^^***^*^***^* *********** 

SSP^’enr J!; ,££'^'■“ 2 '' P£OGrFsSXVE S&HPUTfif! rELL(GRir))WISE 

I 5 l. 5 SiSg 2 i'UJr?,y^|? 4 ;i,« 5 HncjCPR 0 CjPAiPR 0 ^^^^^^^^^^^^^ 

PEAL INKAr( 3 t», 35 ),ri,> 1 ,M D . V.« AGaRA.IU 

real, PATMf,Tt 35 , 3 l 5 


,Ttl 


COMMON/ A! /PXHIM-Py^*x»^ 
coMMor4/A?/L,w,Mjnx^Dy,i 
P£AL HDWB3lH3i,i5) 

XHTECyE^ ^^UXRS/UXC5,PX»^IW,PyHIM 
XMIN,^aM,XMAXry^Uy 
iwTurtfc:p ^i€^pTs 
TNTtiaLR PHAXRS, PWAXCS 
npRN r LIM ITS? 1 , FTbF:=''i->ATA ,DAT’''I 

open (Ur^iT = ?2,PTLiF='c:r.f=:. DAT' 1 
0PEfjCUNUT = 2-i,FILE = 'PAT,DAT^) 


MaKPTS=3^ 

T.=200.?5 
M=0,O0P 
^=•5.028 
THS9. 

read (21 ,*) Wih’Xr ^AXX,MI^’Y,MAXY 
READ(21 f ♦}MAXRS,wiax6s 

READ(2?,») ((XH^AICI,J) ,3=1 , M aXCS ) , 1=1 , kAxR S) 

MAIN PRO^Khm-tl x,"i A ^ 

DX = CMAXX-4TNX)/(MAXCS-1.) 

Dy=C^AXX-HTNy}/CWAXRS-l} 

T«TIT.I7.AmN 

TnVATjID EATRY is denoted by TNVbO 

INVLDs-1 

PhXRS=«AXKS-1 

pmxry=maxc:>-i 

************ t***i^*^***:^******i^^*T^ 

TNTTXDIZIM^. PDWdnUNDRY MATRIX FDR INVALID ENTRIES 
********************************************************** 
DQ 10 T=1,PMAXC5 
DO 10 Jsl.MAXPTS 

ROWBOOCI, J)=Ti>jVbn 
CONTINUE 

INlTLTZTN'; ^ATCH MATRIX 
DU 20 1 = 1 ^MAXpT'S 
DO 20 J=1,MAXPTS 

PAXMATCI,J) = |NVLD'»'X. 

CONXXNUE 

DO 30 I=1,P»1AXRS 
DO 40 t si,PMAXCS 

COf^ING BDJNDRy INFORMATION FOR PATCBCi^Jl 
DD 50 K=t,MAXPTS 
PM'MATtl^K) = R0WR0U(3,K) 

CDNXT4US 

CQPING FRD^ IMPUT MATRIX 
PXMI'J =XMIN+CT-1)»DX 
PYAT>i sYMlN+(J-l)*Dy 
Xl=:(»XNirUDX/2) 

YlsPYMTN 
X2=PXMTM 
y2 = ?YMI.M + 0Y/2 
X3 = X1 
y3=i^2 

X4 = PXmM tDX 

Y4sif2 

X5 = X1 

Y5=pyNTN+DY 

PATMAT (l,i) = TN«AT(T,J) 

PA'f^AT C1.MAXPT5) = ^ I r iJ + U 

PAT'UT (MAXP'TS,!) = TN^ATCTtl ,vl) 

PAT'^AT ( v\AXPTS,MAXPT 5 ) = TNM ATCTtl , J + 1 ) 

PATMAXC i, CMAXPTS+l)/2)-L+N^XltN4yi 

PATMATC(’^AXPTS+1)/2,1 ) = L+.HtX2 + N»Y2 

PATMATC C’^AXPTS+X V2, (MAXPrS^l ) /2)=L + M4fp4-N4' Y3 

PATmATC C«AXPTStl)/2.MAXPTS2=LTM4X4+N»Y4 

PATMATCMAXPlS,CMAXPTwStl)/2)=L + M*X5TN*Y5 
CALL PATPRD C PATMAT, M AXPTS) 

DD 60 TM=i,MAXPTvB 
PD 60 TNslrMAXPTS 
WRIIF C24,t) PATWATCTM^IN) 

cdnitnue 

copy TNG LA$X COLUMN TO FIRST COLUMN 
DO 70 TM =1,VAXPTS 

PATMATf ). = PA3 wATCTM^MAXRTS) 

CONTINUE 

DO SO IMsi ^MAKPXS 
no SO TN~2,MAKFXS 
PATMAT CTVi,lN}=lNVI.l) 


HO rpNTTNNE 

4 n ' rOft.TUitiF 



yo 

S i> 


2i>0 


210 




DQ 90 Tt\ = l,MAKPrS 

paTmat c I /t-o =Tfjvi.,n 

CDVTINOE 

continue: 

Slop 

fnd 

sifBRDUTif'-Js: rn prhcfss p patch 

SUbROUTIt^F PATPP^f PAT^^ aI' HAXpTb) 
DIM£HST0‘I PATMAIf 3^,33) 
integer rj-v-nh ' 

Real l,h.\ 

COMMON/ A Vf.rW,»,,nX,DV,Tri 
DO 43(/P5 lTARZ;rj = i 4 

D £ * M * ! n “ P n W ki n i T T » L> T ; 


X uo 

NERpTSs 0 

IF ?H£k'pTS .FU. CO GO 

X 01 

TFfRDNNu.EO.4) GO TU 

no 

N£WpTS= 0 

DQ 110 Is 1, MAXPTS 
lOsI^l 

Call Rproc ci^^runno, 
continue 

X20 

DO 170 J= 1, MAXPTS 
IQDsJ41, 

CM.b CPRDCdOl). 
continue 

C 

GO TO 100 

Go TO 101 

TyPiJ*- .RUNNO 

continue 

4309*5 

X50 

return 


E^’D 

SiiRRDUTI'^e XD PRHCFSS RCiWS IN PATMAX 

SyRHOUTINE HPRncrRnWNQrRHNND.PATMAT.MAXpTS) 
ni^^tNSTQ'J PATMAr(33,33) 
integer ronno^rowno 
real 

CO«MON/A! / PX*^XN,P^mim 

common/ A? /E,M,N, nX,Dy,TH 

OFFSETS C^AXPISf l5/?. 

DO 200 1= 1, RUNNO 
nFFSFXs rOFFSET^-l) /?. 


CQNXTNME 

no 2J0_J 


MU ijjv j= tj|MAX?TS*“2 

IF (PATMArCftD«N3,j) .EO, INVT.D) GO ID 2l 
IF CPATMAXCROaNO^J + I) .EO. iNVf.D) GOTO 2. 

IF (PATMArUo^NO; J+2) .EQ, TNVLD) 30 Tq 
RFDFIs pa rMATtROi^NO,.!)- PAT.ViAT(R0.^N3, J + i ) 
RFDPPs PAXMATCROWNO^Otl) - P ATMAT C ROWNO 
SDFs kFDFl- RFDF2 
IF (vSDF .UF. XH) GO TO 210 
NEWPTSs NE*^PTS+ ? 

Ji C ents t 
J 2CENTS nCEWX + 1 

RXICFNS PXMINf (ROWND -!)♦ (DX/MAXPTS) 

RX2CEKS RXICFN 

RIlCFNs PX^lNt Cny/MAXPTS)=» (JlCENT -1) 
RY2CFNS RtiZFrA 

TFTPATMATCRO'VNO, JlCENT) .EO, INVLD) GOTq 210 
PATmATCROWNJ, JtCF»<;T)= I.. + M »RX1CEN +n^cryicen) 
PATMAT(RnWNO,02CEivT)s L+MfRX2CFN + N»RY2CEN 

continue 

RETURN 

END 

» f !»; ♦ 

Subroutine xo 

4 t » 

SUPROUTINE CP 
COMMnH/Al/'’XM 

C0MM0N/A2/bt-^ 
dimension pa 
integer RJNiO 
offset sCmAXP 

DO 270 J = Ir 
OFFSEI s (OFE 

cunttrue 

no 230 Xsl,C?* 

TFCPATMATCCOT/ 

TFf PAXMATCCori 
TFfPATMATfCOD 
CDFl sPATHAK 


0 

252. 

“ 210 

J + 7} 


44 4t4)|c444*44444*444*4*t»4»4f4^4» 
PROCFSS COEOUMNS iN PATMAX 
444t4***444t4«4*4.»444^444444 44 44 4# 

ROCCCOLMQ,RDNNO,PAXMAT,MAXpTS) 

IN pyMiN 

rMATC3i,33} 

xs4n/2 

RUNua 

SFX+i)/? 

AXPTS-2) 

Nn,T).ED,TNVLn) GOTO 230 
NO.T + n.EO.TNVLD) GOfO230 
Mn,T + ?).EO.INVL.D) goto 23? 

CObNU^D-PATMATCCObNOrl + l) 


CDF2 sPATMAX(CObNO^ I + ll-PATMArrcObNO^ 14-2) 
rsnp s rr jFi ' . 



if 3D 


TlCE^^T 

T2C£'JT 

CXICFN 

CX7CFN 

CyiCFN 

CY2CFN 


JF (CSDF,bE:tTfi)‘'ho in 230 
NEWPTS s:DeWPTS + 2 
= T 

sIlCSMI+l 

rPX^iM +(C0LNn-l)»(r)X/MRXPr5) 

“ C X 1 C F! * 

= +(nY/MRXPTS)*(IlCE^I-l) 

^- --• ♦ CnX/MaXPTS)*(T2CE^JT-l) 

GDTD 

aI.LsCaICK!’^ 

X2C=CX2CEn 

YICsCYICEm 

Y2C=CY2CFiN 

PATrtAKCnb^D, UCFfn=r.tM*XlCtM>t:yiC 
PATMAKC'^LN'D, l2CFu)=r,.+ «=»:X2C + M*Y2C 

CONTTNDE 4 t 

RETURN 
F i\ D 


230 



u u 


77 

88 

99 

109 

110 


s ‘cn;ir‘L:s:r’*“"^‘-“"^ 

READ (21,*) SS,U,TPS,CSL 
READ (21,*) (X(I),I=i,TPS) 
read (21,*) (Y(I),I=1,tPS) 

DO 77 1=1, TPS 

ZDTM (I) =L+11*X(I)+N*Y(I) 

VIJ =0. 

DO 88 1=1, TPS 
VIJ =VIJ +ZDTM(I) 

ZPDM = VIJ/TPS 
DO 99 1=1, TPS 
ZCOR (I) = ZPDM -ZDTM (I) 

DO 109 I =1,TPS 

AREA (I) = ZCOR(I) *(U+SS*ZCOR(I)) 

DO 110 K= 1, TPS-1 

AREACE+1))«CSL 

URITE (24,*) (VOL(K),K=l, TPS-1) 

STOP 

END 



